#include "HY1C_out.h"
/* =========================================================== */
/* Module loadl1.c                                             */
/*                                                             */
/* Functions to fill a level-1b file with precomputed          */
/* atmospheric and masking data.                               */
/*                                                             */
/* Note: due to filtering capabilities of MSl12, it can not be */
/* assumed that the same l1rec pointer will be passed in later */
/* calls to this function, so all fields must be reloaded. In  */
/* addition, it is possible for MSl12 to process a L1B file    */
/* containing data from different time periods, so earth-sun   */
/* distance can not be assumed to be constant.                 */
/*                                                             */
/* Written By:                                                 */
/*                                                             */
/*     B. A. Franz                                             */
/*     SAIC General Sciences Corp.                             */
/*     NASA/SIMBIOS Project                                    */
/*     February 1999                                           */
/*                                                             */
/* =========================================================== */

#include <stdio.h>
#include "l12_proto.h"
#include "smi_climatology.h"



float *twoD2oneD(float **input,int iscan,int length){
    float *data;
    data=(float *)malloc(sizeof(float)*length);
    int i=0;
    for(i=0;i<length;i++){
        data[i]=input[iscan][i];
    }
    return data;
}


float *threeD2oneD(float ***input,int iscan,int length,int nbands){
    float *data;
    data=(float *)malloc(sizeof(float)*length*nbands);
    int i=0;
    int j=0;
    for(i=0;i<length;i++){
        for(j=0;j<nbands;j++){
            data[i*nbands+j]=input[j][iscan][i];
        }

    }
    return data;
}

int loadl1( filehandle *l1file, instr *input, l1str *l1rec)
{
    static double radeg = RADEG;

    static int   landmask = 0;
    static int   bathmask = 1;
    static int32_t  sensorID = -999;

    static int32_t  *Lambda;
    static int32_t    nwvis;
    static float *Fobar;
    static float *Tau_r;
    static float *k_oz;
    static float *aw;
    static float *bbw;

    int status;
    int32_t ip, ipb, ib, iw, ix;
    double  hour;
    double  esdist;
    int32_t    nbands   = l1rec->nbands;
    int32_t    nbandsir = l1rec->nbandsir;

    double *rvs;
    double temp;

    if (sensorID != l1rec->sensorID) {

        sensorID = l1rec->sensorID;

        if (rdsensorinfo(l1rec->sensorID,l1rec->input->evalmask,NULL,NULL) != l1rec->nbands) {
            HY1C_out("-E- %s : Unable to read sensor file\n",__FILE__);
            exit(1);
        }
        rdsensorinfo(l1rec->sensorID,l1rec->input->evalmask,"Lambda",(void **) &Lambda);
        rdsensorinfo(l1rec->sensorID,l1rec->input->evalmask,"Fobar", (void **) &Fobar);
        rdsensorinfo(l1rec->sensorID,l1rec->input->evalmask,"Tau_r", (void **) &Tau_r);
        rdsensorinfo(l1rec->sensorID,l1rec->input->evalmask,"k_oz",  (void **) &k_oz);
        rdsensorinfo(l1rec->sensorID,l1rec->input->evalmask,"aw",    (void **) &aw);
        rdsensorinfo(l1rec->sensorID,l1rec->input->evalmask,"bbw",   (void **) &bbw);
        nwvis = rdsensorinfo(l1rec->sensorID,l1rec->input->evalmask,"NbandsVIS",  NULL);

        HY1C_out("Loading land mask file from %s\n", input->land);

        if (b128_msk_init(input->land,landmask) != 0) {
            HY1C_out("-E- %s : Unable to initialize land mask\n",__FILE__);
            exit(1);
        }

        HY1C_out("Loading bathymetry mask file from %s\n", input->water);

        if (b128_msk_init(input->water,bathmask) != 0) {
            HY1C_out("-E- %s : Unable to initialize bath mask\n",__FILE__);
            exit(1);
        }

        HY1C_out("Loading ice mask file from %s\n", input->icefile);

        if (ice_mask_init(input->icefile,(int)(*l1rec->year),
                          (int)(*l1rec->day), input->ice_threshold) != 0) {
            HY1C_out("-E- %s : Unable to initialize ice mask\n",__FILE__);
            exit(1);
        }

    }

    /* Get correction for Earth-Sun distance and apply to Fo  */
    esdist = esdist_(l1rec->year,l1rec->day,l1rec->msec);
    l1rec->fsol = pow(1.0/esdist,2);

    for (iw=0; iw<nbands; iw++) {

        l1rec->iwave [iw] = Lambda[iw];
        l1rec->fwave [iw] = (float) Lambda[iw];
        l1rec->Tau_r [iw] = Tau_r [iw];
        l1rec->k_oz  [iw] = k_oz  [iw];
        l1rec->aw    [iw] = aw    [iw];
        l1rec->bbw   [iw] = bbw   [iw];
        l1rec->Fobar [iw] = Fobar [iw];
        l1rec->Fo    [iw] = l1rec->Fobar[iw] * l1rec->fsol;

        get_f0_thuillier_ext(l1rec->iwave[iw],BANDW,&l1rec->Fonom[iw]);
    }

    /* Apply vicarious cross-calibration gains */

    for (ix=0; ix<input->xcal_nwave; ix++) {
        if ((input->xcal_opt[ix] & XCALRVS) != 0) {
            if ((ib = bindex_get(input->xcal_wave[ix])) < 0) {
                HY1C_out("-E- %sline %d: xcal wavelength %f does not match sensor\n",
                     __FILE__,__LINE__,input->xcal_wave[ix]);
                exit(1);
            };
            rvs = get_xcal(l1rec, XRVS, l1rec->iwave[ib]);
            for (ip=0;ip<l1rec->npix; ip++) {
                ipb  = ip*NBANDS+ib;
                if (l1rec->Lt[ipb] > 0.0 && l1rec->Lt[ipb] < 1000.0)
                    l1rec->Lt[ipb] /= rvs[ip];
            }
        }
    }

    for (ip=0; ip<l1rec->npix; ip++) {

        /* Apply vicarious calibration */
        for (iw=0; iw<nbands; iw++) {
            ipb = ip*NBANDS+iw;
            if (l1rec->Lt[ipb] > 0.0 && l1rec->Lt[ipb] < 1000.0) {
                l1rec->Lt[ipb] *= input->gain   [iw];
                l1rec->Lt[ipb] += input->offset [iw];
	    }
        }

        /* Enforce longitude convention */
        if (!l1rec->navfail[ip]) {
            if (l1rec->lon[ip] < -180.)
                l1rec->lon[ip] += 360.0;
            else if (l1rec->lon[ip] > 180.0)
                l1rec->lon[ip] -= 360.0;
        }

        /* Set land/bathymetry flags */
        if (b128_msk_get(l1rec->lat[ip],l1rec->lon[ip],landmask) != 0 && input->format != FMT_L3BIN) {
            l1rec->land[ip] = ON;
        }
        if (!l1rec->land[ip] &&
            b128_msk_get(l1rec->lat[ip],l1rec->lon[ip],bathmask) != 0) {
            l1rec->swater[ip] = ON;
        }
        if (!l1rec->land[ip] &&
            ice_mask(l1rec->lon[ip],l1rec->lat[ip]) != 0) {
            l1rec->ice[ip] = ON;
        }

        /* Force all pixels to water pixels */
        if (input->proc_ocean == 2) {
            l1rec->land[ip] = OFF;
        }

        /* Adjust lon/lat and view geom for terrain height */
        if (!l1file->terrain_corrected) {
          if (input->proc_land && !l1rec->navfail[ip]) {
//            get_dem_height_(input->demfile,
//                            &l1rec->lon [ip],
//                            &l1rec->lat [ip],
//                            &l1rec->senz[ip],
//                            &l1rec->sena[ip],
//                            &l1rec->height[ip],
//                            (int32_t *) &status);
            if (status != 0) {
              HY1C_out("-E- %s line %d: error reading DEM file.\n",
                     __FILE__,__LINE__);
              return(1);
            }
          } else
            l1rec->height[ip] = 0.0;
        }

        /* Set sea surface temperature and salinity, and seawater optical properties */
        for (iw=0; iw<nbands; iw++) {
            ipb = ip*NBANDS+iw;
            l1rec->sw_n [ipb] = 1.334;
 	    l1rec->sw_a [ipb] = l1rec->aw [iw];
            l1rec->sw_bb[ipb] = l1rec->bbw[iw];
	}
        if (input->outband_opt >= 2) {
            for (iw=0; iw<nwvis; iw++) {
                ipb = ip*NBANDS+iw;
 	        l1rec->sw_a [ipb] = aw_spectra (l1rec->fwave[iw],BANDW);
                l1rec->sw_bb[ipb] = bbw_spectra(l1rec->fwave[iw],BANDW);
    	    }
	}
        if (!l1rec->land[ip]) {
            l1rec->sstref[ip] = get_sstref(input->sstfile,l1rec->lon[ip],l1rec->lat[ip],(int)(*l1rec->day));
            l1rec->sssref[ip] = get_sssref(input->sssfile,l1rec->lon[ip],l1rec->lat[ip],(int)(*l1rec->day));
            if (l1rec->sstref[ip] > BAD_FLT && l1rec->sssref[ip] > BAD_FLT && input->seawater_opt > 0) {
              for (iw=0; iw<nwvis; iw++) {
                ipb = ip*NBANDS+iw;
 	        l1rec->sw_n [ipb] = seawater_nsw(l1rec->fwave[iw],l1rec->sstref[ip],l1rec->sssref[ip],NULL);
 	        l1rec->sw_bb[ipb] = seawater_bb (l1rec->fwave[iw],l1rec->sstref[ip],l1rec->sssref[ip]);
	      }
	    }               
        } else {
            l1rec->sstref[ip] = BAD_FLT;
            l1rec->sssref[ip] = BAD_FLT;
	}
        seawater_set(l1rec);

        /* Compute relative azimuth */
        l1rec->delphi[ip] = l1rec->sena[ip] - 180.0 - l1rec->sola[ip];
        if (l1rec->delphi[ip] < -180.)
            l1rec->delphi[ip] += 360.0;
        else if (l1rec->delphi[ip] > 180.0)
            l1rec->delphi[ip] -= 360.0;

        /* Precompute frequently used trig relations */
        l1rec->csolz[ip] = cos(l1rec->solz[ip]/radeg);
        l1rec->csenz[ip] = cos(l1rec->senz[ip]/radeg);

        /* Scattering angle */
        temp   = sqrt((1.0-l1rec->csenz[ip]*l1rec->csenz[ip])*(1.0-l1rec->csolz[ip]*l1rec->csolz[ip]))
               * cos(l1rec->delphi[ip]/radeg);
        l1rec->scattang[ip] = acos(MAX(-l1rec->csenz[ip]*l1rec->csolz[ip] + temp,-1.0))*radeg;
    }

    /* add ancillary data */
    l1rec_HY1C *tmp_hic;
    if ( setanc(l1rec, input,tmp_hic) != 0 )
        return(1);

    if (input->windspeed > -999)
        for (ip=0; ip<l1rec->npix; ip++)
            l1rec->ws[ip] = input->windspeed;
    if (input->windangle > -999)
        for (ip=0; ip<l1rec->npix; ip++)
            l1rec->wd[ip] = input->windangle;
    if (input->pressure > -999)
        for (ip=0; ip<l1rec->npix; ip++)
            l1rec->pr[ip] = input->pressure;
    if (input->ozone > -999)
        for (ip=0; ip<l1rec->npix; ip++)
            l1rec->oz[ip] = input->ozone;
    if (input->watervapor > -999)
        for (ip=0; ip<l1rec->npix; ip++)
            l1rec->wv[ip] = input->watervapor;
    if (input->relhumid > -999)
        for (ip=0; ip<l1rec->npix; ip++)
            l1rec->rh[ip] = input->relhumid;


    /* add atmospheric components that do not depend on Lt */

    for (ip=0; ip<l1rec->npix; ip++) {

        /* ------------------------------------------------ */
        /* Ocean processing                                 */
        /* ------------------------------------------------ */
        if ((input->proc_ocean != 0) && !l1rec->land[ip]) {

  	    if (input->ocrvc_opt == 0)
                atmocor1(l1rec,ip);

            /* set the smile_delta field */
            radcor(l1rec,ip,0);

            /* set polarization correction */
            polcor(l1rec,ip);

            /* add surface reflectance */
            get_rhos(l1rec,ip);

            /* assign uncertainty on Lt if not already set by sensor-specific i/o */
            for (iw=0; iw<nbands; iw++) {
                ipb = ip*NBANDS+iw;
                if (l1rec->Lt_unc[ipb] < BAD_FLT+1) {
	            l1rec->Lt_unc[ipb] = (l1rec->Lt[ipb]-l1rec->Lr[ipb])*input->gain_unc[iw];
                }
            }
        }

        /* ------------------------------------------------ */
        /* Land Processing                                  */
        /* ------------------------------------------------ */
        else if (input->proc_land && l1rec->land[ip]) {
            atmocor1_land(input,l1rec,ip);
            radcor(l1rec,ip,1);
            get_rhos(l1rec,ip);
        }

        /* ------------------------------------------------ */
        /* General Processing                               */
        /* ------------------------------------------------ */
        else {
            for (ib=0; ib<NBANDS; ib++) {
                ipb = ip*NBANDS+ib;
                l1rec->Lr      [ipb] = 0.0;
                l1rec->t_sol   [ipb] = 1.0;
                l1rec->t_sen   [ipb] = 1.0;
                l1rec->tg_sol[ipb] = 1.0;
                l1rec->tg_sen[ipb] = 1.0;
                l1rec->t_o2    [ipb] = 1.0;
                l1rec->t_h2o   [ipb] = 1.0;
                l1rec->polcor  [ipb] = 1.0;
                l1rec->radcor  [ipb] = 0.0;
            }
	}

    }

    /* set masks and flags */
    if ( setflags(input, l1rec) == 0 )
        return(1);

    setflagbits(1,l1rec,NULL,-1);

    return (0);
}


float HY1C_interp(float x1, float   y1, float  x2, float  y2, float  x)
{
    return y1 + (y2 - y1)/(x2 - x1) *(x - x1);
}

void HY1C_spline(float *in, int in_size, float *out, int out_size)
{
    int rate = out_size / in_size;
    int other = out_size % in_size;
    // 内插策略
    // 如果other 为0，整除，不用处理
    // 如果other 不为0，最后一步内插成多个点，确保最值是一致的

    int my_in_size = in_size - 1;
    int my_out_size = (in_size - 1)*rate;
        int last_in_size = 1;
    int last_out_size = rate;
    int i=0,j=0;

    if (other != 0)
    {
        last_out_size = rate + other;
    }

    if (rate == 0)
        return;
    for ( i = 0; i<my_out_size; i++)
    {
        j = i/rate;
        out[i] = HY1C_interp(j, in[j], j + 1, in[j + 1], i/10.0);

    }

    for ( i = my_out_size; i<out_size; i++)
    {
        j = my_in_size;
        out[i] = HY1C_interp(j, in[j-1], j + 1, in[j ], i/10.0);
    }


}



int loadl1_new(filehandle *l1file,instr *input,l1str *l1rec,l1rec_HY1C *data_hic,int iscan)
{
    static double radeg = RADEG;
    static int   landmask = 0;
    static int   bathmask = 1;
    static int32_t  sensorID = -999;
    static int32_t  *Lambda;
    static int32_t    nwvis;
    static float *Fobar;
    static float *Tau_r;
    static float *k_oz;
    static float *aw;
    static float *bbw;
    int status;
    int32_t ip, ipb, ib, iw, ix;
    double  hour;
    double  esdist;
    int32_t    nbands   = 8;
    int32_t    nbandsir = 2;
    double *rvs;
    double temp;

    l1rec->year=year_hy1c_tp ;
    l1rec->day=day_hy1c_tp ;
    l1rec->msec=data_hic->msec;

    if (sensorID != l1rec->sensorID) {

        sensorID = l1rec->sensorID;

        if (rdsensorinfo(l1rec->sensorID,l1rec->input->evalmask,NULL,NULL) != l1rec->nbands) {
            HY1C_out("-E- %s : Unable to read sensor file\n",__FILE__);
            exit(1);
        }
        rdsensorinfo(l1rec->sensorID,l1rec->input->evalmask,"Lambda",(void **) &Lambda);
        rdsensorinfo(l1rec->sensorID,l1rec->input->evalmask,"Fobar", (void **) &Fobar);
        rdsensorinfo(l1rec->sensorID,l1rec->input->evalmask,"Tau_r", (void **) &Tau_r);
        rdsensorinfo(l1rec->sensorID,l1rec->input->evalmask,"k_oz",  (void **) &k_oz);
        rdsensorinfo(l1rec->sensorID,l1rec->input->evalmask,"aw",    (void **) &aw);
        rdsensorinfo(l1rec->sensorID,l1rec->input->evalmask,"bbw",   (void **) &bbw);
        nwvis = rdsensorinfo(l1rec->sensorID,l1rec->input->evalmask,"NbandsVIS",  NULL);

//        HY1C_out("lambda=%d\n",Lambda[2]);
//        HY1C_out("lambda=%d\n",Lambda[3]);
//        HY1C_out("lambda=%d\n",Lambda[4]);
//        HY1C_out("lambda=%d\n",Lambda[5]);
        HY1C_out("Loading land mask file from %s\n", input->land);


        if (b128_msk_init(input->land,landmask) != 0) {
            HY1C_out("-E- %s : Unable to initialize land mask\n",__FILE__);
            exit(1);
        }

        HY1C_out("Loading bathymetry mask file from %s\n", input->water);

        if (b128_msk_init(input->water,bathmask) != 0) {
            HY1C_out("-E- %s : Unable to initialize bath mask\n",__FILE__);
            exit(1);
        }

        HY1C_out("Loading ice mask file from %s\n", input->icefile);

        if (ice_mask_init(input->icefile,(int)(year_hy1c_tp[0]),
                          (int)(day_hy1c_tp[0]), input->ice_threshold) != 0) {
            HY1C_out("-E- %s : Unable to initialize ice mask\n",__FILE__);
            exit(1);
        }

    }

    /* Get correction for Earth-Sun distance and apply to Fo  */

    esdist = esdist_(year_hy1c_tp ,day_hy1c_tp ,(data_hic->msec));
    l1rec->fsol = pow(1.0/esdist,2);

    for (iw=0; iw<nbands; iw++) {

        l1rec->iwave [iw] = Lambda[iw];
        l1rec->fwave [iw] = (float) Lambda[iw];
        l1rec->Tau_r [iw] = Tau_r [iw];
        l1rec->k_oz  [iw] = k_oz  [iw];
        l1rec->aw    [iw] = aw    [iw];
        l1rec->bbw   [iw] = bbw   [iw];
        l1rec->Fobar [iw] = Fobar [iw];
        l1rec->Fo    [iw] = l1rec->Fobar[iw] * l1rec->fsol;

        get_f0_thuillier_ext(l1rec->iwave[iw],BANDW,&l1rec->Fonom[iw]);
    }

    /* Apply vicarious cross-calibration gains */

//    for (ix=0; ix<input->xcal_nwave; ix++) {
//        if ((input->xcal_opt[ix] & XCALRVS) != 0) {
//            if ((ib = bindex_get(input->xcal_wave[ix])) < 0) {
//                HY1C_out("-E- %sline %d: xcal wavelength %f does not match sensor\n",
//                     __FILE__,__LINE__,input->xcal_wave[ix]);
//                exit(1);
//            };
//            rvs = get_xcal(l1rec, XRVS, l1rec->iwave[ib]);
//            for (ip=0;ip<l1rec->npix; ip++) {
//                ipb  = ip*NBANDS+ib;
//                if (l1rec->Lt[ipb] > 0.0 && l1rec->Lt[ipb] < 1000.0)
//                    l1rec->Lt[ipb] /= rvs[ip];
//            }
//        }
//    }

//    l1rec->Lt=threeD2oneD(data_hic->Lt,iscan,data_hic->npix,8);
//    float *tmp_lat=twoD2oneD(data_hic->lat,iscan,data_hic->npix);
//    float *tmp_lon=twoD2oneD(data_hic->lon,iscan,data_hic->npix);

//    float *tmp_solz=twoD2oneD(data_hic->solz,iscan,data_hic->npix);
//    float *tmp_sola=twoD2oneD(data_hic->sola,iscan,data_hic->npix);
//    float *tmp_senz=twoD2oneD(data_hic->senz,iscan,data_hic->npix);
//    float *tmp_sena=twoD2oneD(data_hic->sena,iscan,data_hic->npix);
    float tmp_lat[data_hic->npix2];
    float tmp_lon[data_hic->npix2];
    float tmp_solz[data_hic->npix2];
    float tmp_sola[data_hic->npix2];
    float tmp_senz[data_hic->npix2];
    float tmp_sena[data_hic->npix2];
    // Lt is 0
    l1rec->Lt=(data_hic->Lt+(iscan*8*data_hic->npix));
    // l1rec->Lt=(data_hic->Lt+0*data_hic->nscans*data_hic->npix+(iscan*data_hic->nscans));
    l1rec->lat=(data_hic->lat+(iscan*data_hic->npix));
    l1rec->lon=(data_hic->lon+(iscan*data_hic->npix));
    l1rec->solz=(data_hic->solz+(iscan*data_hic->npix));
    l1rec->sola=(data_hic->sola+(iscan*data_hic->npix));
    l1rec->senz=(data_hic->senz+(iscan*data_hic->npix));
    l1rec->sena=(data_hic->sena+(iscan*data_hic->npix));
    int o=0;
    for(o=0;o<data_hic->npix2;o++){
        tmp_lat[o]=*(l1rec->lat+o);
        tmp_lon[o]=*(l1rec->lon+o);
        tmp_solz[o]=*(l1rec->solz+o);
        tmp_sola[o]=*(l1rec->sola+o);
        tmp_senz[o]=*(l1rec->senz+o);
        tmp_sena[o]=*(l1rec->sena+o);
    }

    float check=tmp_lat[5];

    float tmp_index[data_hic->npix2];
    int ddd=0;
    for(ddd=0;ddd<data_hic->npix2;ddd++){
        tmp_index[ddd]=ddd+1;
    }
    float tmp_index_npix[data_hic->npix];

    for(ddd=0;ddd<data_hic->npix;ddd++){
        tmp_index_npix[ddd]=ddd+1;
    }
//    float check=data_hic->Lt[0][0][data_hic->npix-1];
    if(data_hic->npix!=data_hic->npix2){
        /*lspline(tmp_index,tmp_lat,data_hic->npix2,tmp_index_npix,l1rec->lat,data_hic->npix);
        lspline(tmp_index,tmp_lon,data_hic->npix2,tmp_index_npix,l1rec->lon,data_hic->npix);
        lspline(tmp_index,tmp_solz,data_hic->npix2,tmp_index_npix,l1rec->solz,data_hic->npix);
        lspline(tmp_index,tmp_sola,data_hic->npix2,tmp_index_npix,l1rec->sola,data_hic->npix);
        lspline(tmp_index,tmp_senz,data_hic->npix2,tmp_index_npix,l1rec->senz,data_hic->npix);
        lspline(tmp_index,tmp_sena,data_hic->npix2,tmp_index_npix,l1rec->sena,data_hic->npix);
        */

        HY1C_spline(tmp_lat, data_hic->npix2, l1rec->lat, data_hic->npix);
        HY1C_spline(tmp_lon, data_hic->npix2, l1rec->lon, data_hic->npix);
        HY1C_spline(tmp_solz, data_hic->npix2, l1rec->solz, data_hic->npix);
        HY1C_spline(tmp_sola, data_hic->npix2, l1rec->sola, data_hic->npix);
        HY1C_spline(tmp_senz, data_hic->npix2, l1rec->senz, data_hic->npix);
        HY1C_spline(tmp_sena, data_hic->npix2, l1rec->sena, data_hic->npix);

    }

    for (ip=0; ip<l1rec->npix; ip++) {

        /* Apply vicarious calibration */
//        for (iw=0; iw<nbands; iw++) {
//            ipb = ip*NBANDS+iw;
//            if (l1rec->Lt[ipb] > 0.0 && l1rec->Lt[ipb] < 1000.0) {
//                l1rec->Lt[ipb] *= input->gain   [iw];
//                l1rec->Lt[ipb] += input->offset [iw];
//            }
//        }

        /* Enforce longitude convention */
        if (!l1rec->navfail[ip]) {
            if (l1rec->lon[ip] < -180.)
                l1rec->lon[ip] += 360.0;
            else if (l1rec->lon[ip] > 180.0)
                l1rec->lon[ip] -= 360.0;
        }

        /* Set land/bathymetry flags */
        if (b128_msk_get(l1rec->lat[ip],l1rec->lon[ip],landmask) != 0 && input->format != FMT_L3BIN) {
            l1rec->land[ip] = ON;
        }
        if (!l1rec->land[ip] &&
            b128_msk_get(l1rec->lat[ip],l1rec->lon[ip],bathmask) != 0) {
            l1rec->swater[ip] = ON;
        }
        if (!l1rec->land[ip] &&
            ice_mask(l1rec->lon[ip],l1rec->lat[ip]) != 0) {
            l1rec->ice[ip] = ON;
        }

        /* Force all pixels to water pixels */
        if (input->proc_ocean == 2) {
            l1rec->land[ip] = OFF;
        }

        /* Adjust lon/lat and view geom for terrain height */
        if (!l1file->terrain_corrected) {
          if (input->proc_land && !l1rec->navfail[ip]) {
//            get_dem_height_(input->demfile,
//                            &l1rec->lon [ip],
//                            &l1rec->lat [ip],
//                            &l1rec->senz[ip],
//                            &l1rec->sena[ip],
//                            &l1rec->height[ip],
//                            (int32_t *) &status);
            status=0;
              if (status != 0) {
              HY1C_out("-E- %s line %d: error reading DEM file.\n",
                     __FILE__,__LINE__);
              return(1);
            }
          } else
            l1rec->height[ip] = 0.0;
        }

        /* Set sea surface temperature and salinity, and seawater optical properties */
        for (iw=0; iw<nbands; iw++) {
            ipb = ip*nbands+iw;
            l1rec->sw_n [ipb] = 1.334;
        l1rec->sw_a [ipb] = l1rec->aw [iw];
            l1rec->sw_bb[ipb] = l1rec->bbw[iw];
    }
        if (input->outband_opt >= 2) {
//            for (iw=0; iw<nwvis; iw++) {
           for (iw=0; iw<nbands; iw++) {
                ipb = ip*nbands+iw;
            l1rec->sw_a [ipb] = aw_spectra (l1rec->fwave[iw],BANDW);
                l1rec->sw_bb[ipb] = bbw_spectra(l1rec->fwave[iw],BANDW);
            }
    }
        if (!l1rec->land[ip]) {
            l1rec->sstref[ip] = get_sstref(input->sstfile,l1rec->lon[ip],l1rec->lat[ip], (day_hy1c_tp[0]));
            l1rec->sssref[ip] = get_sssref(input->sssfile,l1rec->lon[ip],l1rec->lat[ip],  (day_hy1c_tp[0] ));
            if (l1rec->sstref[ip] > BAD_FLT && l1rec->sssref[ip] > BAD_FLT && input->seawater_opt > 0) {
              for (iw=0; iw<nwvis; iw++) {
                ipb = ip*nbands+iw;
            l1rec->sw_n [ipb] = seawater_nsw(l1rec->fwave[iw],l1rec->sstref[ip],l1rec->sssref[ip],NULL);
            l1rec->sw_bb[ipb] = seawater_bb (l1rec->fwave[iw],l1rec->sstref[ip],l1rec->sssref[ip]);
          }
        }
        } else {
            l1rec->sstref[ip] = BAD_FLT;
            l1rec->sssref[ip] = BAD_FLT;
    }
        seawater_set(l1rec);

        /* Compute relative azimuth */
        l1rec->delphi[ip] = l1rec->sena[ip] - 180.0 - l1rec->sola[ip];
        if (l1rec->delphi[ip] < -180.)
            l1rec->delphi[ip] += 360.0;
        else if (l1rec->delphi[ip] > 180.0)
            l1rec->delphi[ip] -= 360.0;

        /* Precompute frequently used trig relations */
        l1rec->csolz[ip] = cos(l1rec->solz[ip]/radeg);
        l1rec->csenz[ip] = cos(l1rec->senz[ip]/radeg);

        /* Scattering angle */
        temp   = sqrt((1.0-l1rec->csenz[ip]*l1rec->csenz[ip])*(1.0-l1rec->csolz[ip]*l1rec->csolz[ip]))
               * cos(l1rec->delphi[ip]/radeg);
        l1rec->scattang[ip] = acos(MAX(-l1rec->csenz[ip]*l1rec->csolz[ip] + temp,-1.0))*radeg;
    }

    /* add ancillary data */

    if ( setanc(l1rec, input,data_hic) != 0 )
        return(1);

    if (input->windspeed > -999)
        for (ip=0; ip<l1rec->npix; ip++)
            l1rec->ws[ip] = input->windspeed;
    if (input->windangle > -999)
        for (ip=0; ip<l1rec->npix; ip++)
            l1rec->wd[ip] = input->windangle;
    if (input->pressure > -999)
        for (ip=0; ip<l1rec->npix; ip++)
            l1rec->pr[ip] = input->pressure;
    if (input->ozone > -999)
        for (ip=0; ip<l1rec->npix; ip++)
            l1rec->oz[ip] = input->ozone;
    if (input->watervapor > -999)
        for (ip=0; ip<l1rec->npix; ip++)
            l1rec->wv[ip] = input->watervapor;
    if (input->relhumid > -999)
        for (ip=0; ip<l1rec->npix; ip++)
            l1rec->rh[ip] = input->relhumid;


    /* add atmospheric components that do not depend on Lt */
    int hic_npix=data_hic->npix;
    int back_epix=l1rec->epix;
    int back_escan=l1rec->escan;
    l1rec->epix=data_hic->npix-1;
    l1rec->escan=data_hic->nscans-1;
    for (ip=0; ip<hic_npix; ip++) {
//        HY1C_out("ip=%d\n",ip);

        /* ------------------------------------------------ */
        /* Ocean processing                                 */
        /* ------------------------------------------------ */
        if ((input->proc_ocean != 0) && !l1rec->land[ip]) {

        if (input->ocrvc_opt == 0){
            atmocor1(l1rec,ip);
//            HY1C_out("atmocor1\n");
        }


            /* set the smile_delta field */
            radcor(l1rec,ip,0);

            /* set polarization correction */
            polcor(l1rec,ip);

            /* add surface reflectance */
            get_rhos(l1rec,ip);

            /* assign uncertainty on Lt if not already set by sensor-specific i/o */
            for (iw=0; iw<nbands; iw++) {
                // ipb = ip*NBANDS+iw;
                ipb = ip*nbands+iw;
                if (l1rec->Lt_unc[ipb] < BAD_FLT+1) {
                l1rec->Lt_unc[ipb] = (l1rec->Lt[ipb]-l1rec->Lr[ipb])*input->gain_unc[iw];
                }
            }
        }

        /* ------------------------------------------------ */
        /* Land Processing                                  */
        /* ------------------------------------------------ */
        else if (input->proc_land && l1rec->land[ip]) {
            atmocor1_land(input,l1rec,ip);
            radcor(l1rec,ip,1);
            get_rhos(l1rec,ip);
        }

        /* ------------------------------------------------ */
        /* General Processing                               */
        /* ------------------------------------------------ */
        else {
            for (ib=0; ib<nbands; ib++) {
                ipb = ip*nbands+ib;
                l1rec->Lr      [ipb] = 0.0;
                l1rec->t_sol   [ipb] = 1.0;
                l1rec->t_sen   [ipb] = 1.0;
                l1rec->tg_sol[ipb] = 1.0;
                l1rec->tg_sen[ipb] = 1.0;
                l1rec->t_o2    [ipb] = 1.0;
                l1rec->t_h2o   [ipb] = 1.0;
                l1rec->polcor  [ipb] = 1.0;
                l1rec->radcor  [ipb] = 0.0;
            }
    }

    }

    /* set masks and flags */
    if ( setflags(input, l1rec) == 0 )
        return(1);

    setflagbits(1,l1rec,NULL,-1);

    return (0);
}

